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Expansion dynamics of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state 

J. Kajala,^ F. Massel,^ and P. Torma^ 

^Department of Applied Physics, Aalto University School of Science, P.O.Box 15100, FI-00076 Aalto, FINLAND 
^Low Temperature Laboratory, Aalto University School of Science, P.O. Box 15100, FL00076 Aalto, FINLAND 

We consider a two- component Fermi gas in presence of spin imbalance, modelling the system in 
terms of a one-dimensional (ID) attractive Hubbard Hamiltonian initially in presence of a confining 
trap potential. With the aid of the time-evolving block decimation (TEBD) method, we investigate 
the dynamics of the initial state when the trap is switched off. We show that the dynamics of a 
gas initially in the FFLO state is decomposed into the independent expansion of two fluids, namely 
the pairs and the unpaired particles. In particular, the expansion velocity of the unpaired cloud is 
shown to be directly related to the FFLO momentum. This provides an unambiguous signature of 
the FFLO state in a remarkably simple way. 



Ultracold gases have provided experimental verifi- 
cation for several fundamental concepts of quantum 
physics. The direct access to momentum distribution 
and correlations via time-of-flight expansion has played 
a key role in many landmark experiments. Mapping of 
momentum to position after time-of-flight revealed Bose- 
Einstein condensation [1]. Coherence manifesting after 
expansion gave evidence of the phase of the condensate 
[2] , and of the superfluid - Mott insulator transition [3j . 
Inversion of the aspect ratio of an expanding Fermi gas 
revealed hydrodynamic behaviour [j, l5|. A major goal is 
to observe the elusive Fulde-Ferrel-Larkin-Ovchinnikov 
(FFLO) state which is at the heart of understanding 
the interplay between superconductivity and magnetism. 
While experiments in solid state systems [6-8] and ultra- 
cold gases [9] are consistent with the state, an unambigu- 
ous observation is lacking. 

The FFLO phase is characterized by the formation 
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of pairs with nonzero overall momentum p^, [Ijj . This 
property manifests itself through the appearance of an 
oscillating (superconducting) order parameter ^(r) ex 
Aq exp [iqr] , where q is the so-called FFLO momentum 
{h = 1) [12]. It is given by the mismatch between the 
Fermi momenta of N^ spin up and N^ spin down parti- 
cles: q = kp^ — kFi' In contrast, zero-momentum {q = 0) 
pairs give the conventional BCS superconductor physics. 
Despite the lack of genuine long range order, the FFLO 
state has been theoretically predicted to be especially 
stable in one-dimensional systems [13|4l5|, as well as in 
quasi- ID [16|4l8|. Inspired by the first experiments in 
imbalanced atomic Fermi gases [19, 20], various methods 
for detecting the FFLO state in ultracold gases have been 
proposed (see [21] and references therein). In this letter 
we show via exact simulations that time-of-flight expan- 
sion provides a smoking gun signature of the FFLO state 
in one dimension, see Figure [H Expansion and conse- 
quent imaging of densities is a widely used basic tech- 
nique in experiments with ultracold gases. 

The characteristic parameters can be chosen so that 
the lattice model employed is a good approximation for 
a continuum model of a spin-imbalanced gas in ID in 
the strong interaction limit [l6|, |22|, |23| . With appropri- 



FIG. 1: (Color online) Expansion of a spin-density imbalanced 
Fermi gas in the ID FFLO state after the confining potential 
has been switched off. Our exact simulations show an effective 
two- fluid behaviour: pairs and unpaired particles expand with 
different velocities. Remarkably, the expansion velocity of the 
unpaired particles is directly related to the FFLO momentum 
and provides a straightforward way for observing the FFLO 
state. 



ate mapping of the parameters, in the strong interaction 
limit, our predictions are thus relevant for experiments 
in ID potentials as in [9[. 

We consider the ID Fermi-Hubbard Hamiltonian in 
presence of an overall confining potential 



H = U^hi^hii-J ^ 4^Q+i^ + /i.c.+ 



ia='\-,l 



E 



V^fiu 



(1) 

where c\^ {cia) creates (annihilates) a spin a particle at 
the lattice site i G {1, i^}, J is the hopping constant, U is 
the interaction strength between the two species, and Vi 
the strength of the trapping potential. For a harmonic 
potential Vi = Vho{i — C)^, where C denotes the trap 
center, and for a box potential Vi is specified below. 

To obtain the ground state and time evolution of 
the system, we have employed the (essentially) ex- 
act time-evolving block decimation (TEBD) algorithm 
|24| with Schmidt number F = 150 and simulation 
timestep At = 0.02 j. As results from TEBD we obtain 



the single particle densities n^^(t) = (^(t)|cLQ^|^(t)), 

ni^{t) = {^{t)\cl^Ci^\^{t)), the doublon density n^;(t) = 

(<l>(t)|cLcj|Q^Q||$(t)), and the ground state pair corre- 
lation Cij and its momentum transform n/^, which are 
given by 



C^J = (^l4t4;CjiCjtl^). 



n. = ^t^'^'-^''C^: 



(2) 



(3) 



where (<^||<^) describes the quantum mechanical average 
over the state $, i and j are lattice site indices, t is the 
imaginary unit and L is the size of the lattice. 

We have experimented within TEBD paramater ranges 
iV^ = - 40, iV^ = - 40, P = ^^ = 0.024 - 1, 
L = 80 - 320, Vho = 0.02 J - 0.0001 J, and T = 80 - 200, 
where P is the polarization and Vho is the harmonic trap- 
ping strength. Many simulations were performed also us- 
ing a box potential. Single runs with the larger parame- 
ters have been done in order to check that the qualitative 
description of the dynamics stays the same when chang- 
ing the parameters. Indeed, the important characteristics 
of the dynamics are the same in all of the scenarios. 

For simplicity we first consider here the box trap. In 
Figure [2] we show results for 10 up-particles and 6 down- 
particles, with U = —10 J. The ground state is a ID 
FFLO state, characterized by a peak in the pair momen- 
tum distribution n^ that coincides with the definition of 
the FFLO momentum q = kp^ — kpi- Figure [2k shows 
the ground state density profile in which the small oscil- 
lations characterize the FFLO state. However, such del- 
icate features are hard to resolve in ultracold gas experi- 
ments, and therefore other signatures are needed. Figure 
[lla displays the pair momentum correlation function. 

The dynamics after releasing the particles from the po- 
tential show a striking two-fiuid behaviour. Pairs (dou- 
blons) and excess unpaired majority particles expand as 
effectively non- inter acting fiuids, as can be seen from Fig- 
ure [31 Indeed, we have compared the dynamics with 
strong interaction between the spins to the dynamics of 
non-interacting particles and verified that the doublons 
and unpaired particles in the strongly interacting limit 
expand qualitatively just like non-interacting particles 
would (see supplementary information Figures 2-5 for the 
C/ = results [25]). The important difference is, however, 
that the velocity of the expansion is changed with respect 
to the non-interacting case. The doublons expand with 
velocities up to ^sin(/ci?|). And what is crucial for 
our proposal for observing the FFLO state: the unpaired 
particles expand with velocities up to 2Jsin(g), where q 
is the FFLO momentum. In the U ^ J case shown here 
the unpaired particle velocity is larger than the pair one, 
but in general the relative velocities of pairs vs. unpaired 
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FIG. 2: (Color online) a: The density profiles of up (n^) and 
down spins (n^) in the ground state when N^ = 10, Ni = 6, 
U — —10 J, and there is to a good approximation an infinitely 
strong repulsive potential everywhere except the at the lattice 
sites 66-85. b: The pair momentum correlation function Uk 
for the same state. There are peaks at the FFLO momenta 
q = ±(/cFt - kpi) = 0.27r/L. 



particles is not essential for the method since the up and 
down components can be imaged separately, see e.g. [9|. 
A wavefront corresponding to q clearly separates from 
the rest of the unpaired particle cloud during initial dy- 
namics {t = 0— 10 j), after which it moves with a constant 
velocity at the edge of the cloud, see FiglSt,. Therefore, 
by measuring the expansion velocity of the cloud edge 
{vexp)^ e.g. from the maximum gradient of the density 
(see supplementary information Figure 1 [25|), one ob- 
tains the FFLO momentum from 



arcsm 






(4) 



The momentum q obtained in this way from our simula- 
tions is compared in Figure [4^ to the FFLO momentum 
as given by the definition q = kp,^ — kp^i (and by the 
peak in the momentum pair correlation function of the 
ground state, c.f. Figure [SJa). For reference, we show the 
expansion velocities in the cases of a non-interacting gas, 
and a non-FFLO state without pair coherence (discussed 
later in the text). Only in the case of the FFLO state 
does the expansion velocity depend on the imbalance. 
The q extracted from the simulations via Eq.(|l]) matches 
excellently the expected FFLO momentum. 

In the experimentally relevant case of a harmonic trap, 
q can be determined in the same way as discussed above. 
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FIG. 3: (Color online) a: The time development of the dou- 
blon density n^^, corresponding to the ground state shown in 
Figure [21 b: The time development of the unpaired particle 
density n^ — n^^ . 



Figure |4]b shows the comparison of q obtained from the 
edge expansion velocity to the FFLO momentum in a 
harmonic trap. The time evolution is shown in Figure [5l 
In q = kp^ — kpi^ the Fermi momenta are the highest 
momenta in the harmonic oscillator eigenstates of the 
quantum numbers Uf^ and npi^ respectively (see sup- 
plementary information [25]). In the limit of a shallow 
trap and large particle number, one approaches the box 
potential case and a local density approximation (LDA) 
argument can be used for the extraction of kp a !26'] . 

The interpretation of our numerical results in terms 
of the expansion of two fluids is supported by a rigor- 
ous Bethe-ansatz analysis of the one-dimensional Fermi- 
Hubbard Hamiltonian 22|, |23|, l27| . In the strong-coupling 
limit, the system can be described as two weakly inter- 
acting spinless Fermi gases, corresponding to pairs and 
unpaired particles whose maximum group velocities are 
the respective Fermi velocities. We thus expect the time 
evolution of the system to be approximately generated by 
the free-particle Hamiltonians for unpaired particles and 
pairs. In the case of a lattice model, the group velocity for 
each momentum component is Vk = 2 J sin [/c], where the 
constant J depends on the nature of the particle consid- 
ered. In our pair/unpaired particle 2-fluid model stem- 
ming from the Hubbard Hamiltonian, we have J = J 
for the unpaired particles and J = 2J^//7 for the pairs. 
For k < 7r/2, the relation Vk — 2 J sin [/c] thus allows to 
establish a connection between the maximal expansion 
velocity and the Fermi momentum of each component 
(note that for k > ix 12 the maximal expansion velocity is 
given by 2 J). 

Based on our numerical findings, we now assume that 
the unpaired-particle Fermi momentum and the FFLO 
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FIG. 4: (Color online) a: The FFLO momentum q determined 
from the edge expansion velocity of the unpaired cloud, com- 
pared to the q of the ground state, as a function of N^ describ- 
ing the imbalance (N^ — 10). We also show the expansion mo- 
mentum k ^ q obtained in the case of a non-interacting gas, 
and a non-FFLO state without pair coherence, b: Same as a, 
but having initially a shallow harmonic trap {Vho = 0.0003, 
C = 75.5). 
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FIG. 5: (Color online) a: The time development of the dou- 
blon density n^^ with the same parameters as in Figure [2] but 
the initial trap being harmonic [Vho — 0.0003, C — 75.5) in- 
stead of a box trap, b: The time development of the unpaired 
particle density n^ — n^^ . 



momentum share the same value q = kp^ — k^^. This 



is supported by Bethe-ansatz analysis in the continuum 
case [28| and is also intuitive: in the limit of strong in- 
teractions the pair kinetic energies are small and it is 
energetically favourable that the ground state structure 
allows the lowest momentum states up to q to be mainly 
occupied by unpaired majority particles. Thus measuring 
the maximal expansion velocity of the unpaired particles 
gives access to the value of the FFLO momentum. As 
seen in Figure IH this scenario agrees with the numerical 
results, providing a clear signature of the FFLO state. 
The lattice reproduces the continuum-case dynamics for 



low densities [16|, |22| (we have tested also the low density 
limit), and our proposal is thus expected to be suitable 
for the detection of the FFLO state in experiments of the 
type [9[. Moreover, the proposed method should be ro- 
bust with respect to averaging over an array of ID tubes, 
c.f. p. 

Our simulations describe the zero temperature case. 
Quantum Monte Carlo calculations suggest that there is 
a finite-temperature precursor for the FFLO state which 
involves pairing but no FFLO-type correlations [29|. Is 
the signature we propose for the FFLO state distinguish- 
able from any traces of a state with pairing but no FFLO 
correlations? Simulating exact dynamics at finite tem- 
perature is a formidable task and well beyond the state- 
of-the-art. However, to consider pairing without FFLO 
correlations, we have simulated the dynamics of various 
initial states which have the same average densities of un- 
paired and paired particles as the corresponding FFLO 
states, but which possess no correlations between the par- 
ticles or between the pairs. Such an uncorrelated state 
can be written as an ensemble of product states, in which 
each of the product states have completely spatially lo- 
calized pairs and single majority particles. Choosing the 
constituents of the ensemble randomly and giving them 
equal weight, we obtain the results shown in Figure HI 
the expansion velocity does not depend on the polariza- 
tion. It is simply given by 2 J and AJ'^ /U for unpaired 
particles and pairs, respectively. Therefore, observing a 
change of the unpaired particle expansion velocity with 
the polarization, such that it has the functional depen- 
dence g = /cFt ~ ^Fi, is a genuine signature of the FFLO 
correlations, and cannot be achieved for a non-correlated, 
paired state with the same imbalance. When lowering 
the temperature, emergence of a g-dependent expansion 
front on top of a thermal state background would reveal 
that the FFLO state has been reached. 

To conclude, we have shown that like in many classic 
ultracold gas experiments, also in the case of the long- 
sought-for FFLO state, the expansion of the cloud gives 
an exceedingly simple way of determining the nature of 
the initial state. Our exact quantum many-body simula- 
tions of the ID imbalanced gas present a clear two-fluid 
behaviour with characteristic expansion velocities. We 
show that the matching of the expansion velocity of the 
unpaired majority cloud with the expected FFLO mo- 



mentum provides an unambiguous signature of the FFLO 
state. 
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UNCORRELATED STATE DYNAMICS 



SUPPLEMENTARY MATERIAL 

DYNAMICS OF THE GAUDIN-YANG MODEL 
IN THE STRONGLY ATTRACTIVE REGIME 

As shown in [1], the Bethe-ansatz solution for the 
Gaudin-Yang model in the strongly interacting regime 
can be written in terms of the quasi-momenta ki {n for 
unpaired particles, p for pairs) and the spin roots Ap 



^ TT 5n 

kn = T^^ ~r 
L cL 

TT ep 



'-P 



n={l...N^-Ni,} 
p= {1/2,1,3/2... NJ2} 



kp± =^p± 2^c 



(5) 
(6) 

(7) 



where L represents the size of the system and 2c the 
contact interaction strength, while 6 and e are constants 
related to N^/L and N^/L. Note that c has been rescaled 
so that the physical value of the interaction strength is 
given by c = Ji?c/2m. Eqs. (J5]|7|) allow us to rewrite the 
dispersion relation for the system 



E 



in the following form 




E 



2m 



E(2A^cV2)+E^^ 



(8) 



(9) 



From Eq. (|9|), we recognize the pair contribution to the 
kinetic energy (2A^), along with the pair-binding energy 
(c^/2) and the single-particle kinetic energy (/c^). We can 
thus state that since, in the strongly interacting limit, the 
binding energy depends on the interaction energy c only, 
it is independent of the Fermi energy. Remembering the 
definition of Ap, for large interaction (i.e. 1/c ^ ), it 
is possible to rewrite the pair kinetic energy term as 



h 

m V2L' 



^skin = — I :7-r5 1 



(10) 



with s = {1,2...A^|}. The term j^s can be identified 
with the quasi momentum Ks of a free particle in a system 
of size L. If we calculate the velocity of such object 



dE, 



Vs = 



kin 



dhtvs 



flKs 

2m 



(11) 



we can conclude that, in the strong-coupling limit, along 
with N^ — Ni excitations with group velocity hkn/m 
(particles), N^ excitations moving with group velocity 
hi<is/2m. are present. We identify these excitations with 
pairs with binding energy c^/2 and mass 2m. 



The uncorrelated state dynamics originates from the 
dynamics of localized particles released into the lattice, 
similiar to the case previously studied in [2|. The ex- 
pansion wavefronts which emerge have velocities 2 J and 
^jj^ , which is in accordance with our understanding that 
we see the maximum velocities of unpaired particles and 
pairs in the dynamics. Note that localization to a lat- 
tice site corresponds in the Fourier space to employ- 
ing all momenta in the band. The maximum expan- 
sion velocities observed correspond to momenta k = ^ 
since, due to the lattice dispersion E = — 2Jcos(/c), 
Vexp = "^ = 2Jsin(A:), where we have J = J for un- 



paired particles and J = -jj- for pairs. 



DETERMINING THE EDGE EXPANSION 
VELOCITY 

The edge expansion velocity has been determined us- 
ing the numerically obtained density profiles at differ- 
ent times after the release from the trapping potential. 
We are interested in the maximum group velocity of the 
unpaired component, as this gives the FFLO momen- 
tum q. In the simulations it is seen that this unpaired 
wavefront separates from the rest after initial dynamics 
{t = 1 — 10-j) after which it travels at constant velocity at 
the edge of the cloud. The velocity of the edge wavefront 
has been determined by calculating the change of the po- 
sition of the maximum gradient of the unpaired majority 
particle density at the edge. Figure [6] illustrates how this 
has been done. 



COMPARING NON-INTERACTING AND 

STRONGLY INTERACTING EXPANSION 

DYNAMICS 

Figures [7] - [JO] illustrate how the expansion profiles of 
X doublons and Y unpaired particles look qualitatively 
like the expansion profiles of X and Y non-interacting 
particles, respectively, but have different velocities com- 
pared to the non-interacting case. The velocity in the 
strongly interacting regime is 2Jsm{q) for the unpaired 
particles and ^^sin{kFi) for the pairs (down is the mi- 
nority species). In the noninteracting case, the expansion 
velocity of species a is given by 2Jsin(A:Fcr)- 



ADDITIONAL DATA 

Figure [TT] depicts the ground state of the ID FFLO 
superfiuid in a harmonic trap as described in the main 
text. The density is low enough so that the lattice re- 
sult corresponds to the continuum case in the strongly 



interacting regime. Corresponding to this ground state, 
Figure [12] shows the momenta involved in FFLO pair- 
ing in the case of a harmonic trap with smah particle 
number. 

Looking at Figure [121 the q determined from the ex- 
pansion velocity matches kp^ — kri which has been cal- 
culated from noninteracting quantum harmonic oscillator 
eigenstates. To elaborate, for example kp^ for five parti- 
cles has been obtained from the peak at maximum mo- 
mentum in the momentum distribution of the 5th quan- 
tum harmonic oscillator eigenstate (i.e. n = 4, since the 
harmonic oscillator quantum numbers start from n = 0). 
Moreover, the pair momentum correlator Uk matches 
A:(7v-i^_7V^)^ — koi calculated again using the noninteract- 
ing harmonic oscillator eigenstates {koi = for the 1st 
eigenstate which is a Gaussian centered at zero). The 
quantity k(^N^_N^^^ — ko^ corresponds to pairing between 
a down particle at the lowest level (n = 0) and an up 



particle at the level n = N^ — N^^ whereas kp^ — kpi cor- 
responds to a pairing between a down particle at np^ and 
an up particle at np^- For larger particle numbers the 
peak of Uk and kp^ — kp^ (and thus q as obtained from 
the expansion velocity) converge into the same value, as 
verified by our numerical simulatons with larger parti- 
cle numbers (see Figure [14]). In comparison, for the box 
potential, the peak of n/c, kp^ — kp^ and q as obtained 
from the expansion velocity are ahead the same for small 
particle numbers due to the box eigenstates being mo- 
mentum eigenstates (see Figure [T3|) . For larger particle 
numbers, the box potential becomes a reasonable approx- 
imation for the harmonic potential, and thus the values 
of the peak of n/c, kp^ — kp^ and q as obtained from 
the expansion velocity become the same. This explains 
why the local density approximation works also in a har- 
monic trap for determining the FFLO momentum using 
kp^ — kpi = n^ — n^ for large particle numbers. 
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FIG. 6: a: Obtaining the expansion velocity from the position of the maximum gradient of unpaired particle density at the edge. 

13J 

The FFLO wavevector is obtained from q = arcsin(^^|^). Substituting gives q = arcsin(-^)-^ = 0.71 -jj. Particle numbers here 

are Nf = 10, A'4. = 6, and initially the gas was confined in a box potential, b: The same as a but for A^^ = 10, A'4. = 4. The 

16J 
FFLO wavevector in this case is q = arcsin(-^).ij = 0.93-^. 
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FIG. 7: a: The density profile of unpaired particles n^{t) — 
n^l (t) after the release of the gas from the box potential when 
N^ = 10, Ni = S, U = -10.0 J. b: The density profile n^{t) 
of two noninteracting particles released from the box potential 
(A^.^ = 2, Ni — 0). Notice that the time axis is different in the 
two plots, showing how the expansion is qualitatively similiar 
but occurs at a different velocity. 



FIG. 9: a: The density profile n^(t) — n^iit) for parameters 
N^ = 10, N^ = 2, and U = -10.0 J. b: The density profile 

n^{i) for AT-^ = 8 and N^ = 0. 




20 40 60 80 100 120 140 

Lattice site 



I 



I 




100 120 140 



10.3 

0.25 
-0.2 

0.15 
-0.1 

10.05 





I 



0.18 
0.16 
0.14 
0.12 




0.06 
0.04 
0.02 




20 40 60 80 100 120 140 

Lattice site 



I 



FIG. 8: a: The density profile of paired particles n^i{t) after 
the release of the gas from the box potential when N^ = 10, 
Ni = 8, and U = -10.0 J. b: The density profile n^{t) of 
eight noninteracting particles released from the box potential 
(^N^ — 8, Ni — 0). Notice again the scaling of the time axis. 



FIG. 10: a: The density profile n^^{t) for N^ = 10, N^ = 2, 
and U = -10.0 J. b: The density profile n^{t) for N^ ^ 2 and 
iVi = 0. 




FIG. 11: The up and down particle densities n^ and n^. for the 
ground state with harmonic trapping and parameters N^ = 
10, AT^ = 6, L = 150, U = -lO.OJ, and Vho = 0.0003J. 



_^_ Velocity of n^ - n^^, U = -10.0 _ 
_,^ Velocity of n^ - n^^, U = 0.0 




FIG. 12: The FFLO momentum ^ in a harmonic trap as de- 
termined from the edge expansion velocity, the momentum of 
the edge expansion for the noninteracting case, /ci?t ~ kpi as 
obtained from the peaks in momentum distribution of high- 
est occupied (noninteracting) harmonic oscillator eigenstates 
(see text for explanation), /c(iv^_iVi )t5 position of the peaks 
in the pair momentum correlation Uk, and the momentum of 
the edge expansion given by an uncorrelated state. 
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FIG. 14: a: The ground state in a harmonic trap with the 
parameters L = 320, N^ = 40, N^ = 20, Vho = 0.0003, 
and U — —10.0. b: The corresponding time evolution of 
unpaired particles n^{t) — n^i{t) and c: doublons n^i{t). d: 
Pair correlation function rik for the ground state. The peak 
in Uk and q as determined from the edge expansion velocity 
of unpaired particles give the same value q = 0.67;^. 



FIG. 13: The FFLO momentum ^ in a box potential as de- 
termined from the edge expansion velocity, kp^ — kpi as ob- 
tained from the highest occupied noninteracting box eigen- 
states (kpa — ^^^)^ and the position of the peaks in the pair 
momentum correlation Uk- 
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